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Abstract 

With the aim of determining the statistical properties of relativistic turbulence and unveiling novel and non- 
classical features, we present the results of direct numerical simulations of driven turbulence in an ultrarelativis- 
tic hot plasma using high-order numerical schemes. We study the statistical properties of flows with average 
Mach number ranging from ~ 0.4 to ~ 1.7 and with average Lorentz factors up to ~ 1.7. We find that flow 
quantities, such as the energy density or the local Lorentz factor, show large spatial variance even in the sub- 
sonic case as compressibility is enhanced by relativistic effects. The velocity field is highly intermittent, but its 
power-spectrum is found to be in good agreement with the predictions of the classical theory of Kolmogorov. 
Overall, our results indicate that relativistic effects are able to significantly enhance the intermittency of the 
flow and affect the high-order statistics of the velocity field, while leaving unchanged the low-order statistics, 
which instead appear to be universal and in good agreement with the classical Kolmogorov theory. To the best 
of our knowledge, these are the most accurate simulations of driven relativistic turbulence to date. 
Subject headings: Relativistic Hydrodynamics, Turbulence 



1. INTRODUCTION. 

Turbulence is an ubiquitous phenomenon in nature as it 
plays a fundamental role in shaping the dynamics of sys- 
tems ranging from the mixture of air and oil in a car en- 
gine, up to the rarefied hot plasma composing the intergalactic 
medium. Relativistic hydrodynamics is a fundamental ingre- 
dient in the modeling of a number of systems characterized by 
high Lorentz-factor flows, strong gravity or relativistic tem- 
peratures. Examples include the early Universe, relativistic 
jets, gamma-ray-bursts (GRBs), rel ativistic he avy-ion colli- 
sions and core-collapse supernovae (Font 2008). 

Despite the importance of relativistic hydrodynamics and 
the reasonable expectation that turbulence is likely to play 
an important role in many of the systems mentioned above, 
extremely little is known about turbulence in a relativistic 
regime. For this reason, the study of relativistic turbulence 
may be of fundamental importance to develop a quantitative 
description of many astrophysical systems. To this aim, we 
have performed a series of high-order direct numerical simu- 
lations of driven relativistic turbulence of a hot plasma. 

2. MODEL AND METHOD. 

We consider an idealized model of an ultrarelativistic fluid 
with four- velocity = W(l,v l ), where W = (1 — 
ViV 1 )^ 1 / 2 is the Lorentz factor and v l is the three-velocity 
in units where c = 1. The fluid is modeled as perfect and 
described by the stress-energy tensor 

T liV = {p+p)u IJl u l ,+pg l j tv , (1) 

where p is the (local-rest-frame) energy density, p is the pres- 
sure, the four-velocity, and <? M „ is the spacetime metric, 
which we take to be the Minkowski one. We evolve the 
equations describing conservation of energy and momentum 
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in the presence of an externally imposed Minkowskian force 
F^, i.e. V ' V T^ V = F^, where the forcing term is written as 
= F(0, f l ). More specifically, the spatial part of the 
force, /*, is a zero-average, solenoidal, random, vector field 
with a spectral distribution which has compact support in the 
low wavenumber part of the Fourier spectrum. Moreover, f % , 
is kept fixed during the evolution and it is the same for all the 
models, while F is either a constant or a simple function of 
time (see below for details). 

The time component of the forcing term, F°, is zero, so 
that the driving force is able to accelerate fluid elements with- 
out changing their total energy (in the Eulerian frame). On 
the other hand, we impose a minimum value for the energy 
density in the local-rest-frame, p m ; n . This choice is motivated 
essentially by numerical reasons (the very large Lorentz fac- 
tor produced can lead to unphysical point-wise values of p) 
and has the effect of slowly heating up the fluid. Furthermore, 
this floor does not affect the momentum of the fluid and only 
the temperature is increased. From a physical point of view, 
our approach mimics the fact that in the low-density regions, 
the constituents of the plasma are easily accelerated to very 
high Lorentz factors, hence emitting bremsstrahlung radiation 
heating up the surrounding regions. The net effect is that en- 
ergy is subtracted from the driving force and converted into 
thermal energy of the fluid, heating it up. In general p m ; n is 
chosen to be two orders of magnitude smaller than the initial 
energy density, but we have verified that the results presented 
here are insensitive to the specific value chosen for p m - ln by 
performing simulations where the floor value is changed by 
up to two orders of magnitude without significant differences. 

The set of relativistic-hydrodynamic equations is closed by 
the equation of state (EOS) p = hp, thus modelling a hot, 
optically-thick, radiation-pressure dominated plasma, such as 
the electron-positron plasma in a GRB fireball or the matter in 
the radiation-dominated era of the early Universe. The EOS 





Figure 1. Left panel: average Lorentz factor as a function of time for the different models considered. Note that a quasi-stationary state is reached before t ~ 10 
for all values of the driving force. Right panel: logarithm of the Lorentz factor on the (y, z) plane at the final time of model D. Note the large spatial variations 
of the Lorentz factor with front-like structures. The time-averaged PDFs are shown in the lower left corner for the different models considered. 



used can be thought as the relativistic equivalent of the clas- 
sical isothermal EOS in that the sound speed is a constant, 
i.e. <? s = 1/3. At the same time, an ultrarelativistic fluid is 
fundamentally different from a classical isothermal fluid. For 
instance, its "inertia" is entirely determined by the tempera- 
ture and the notion of rest-mass density is lost since the latter 
is minute (or zero for a pure photon gas) when compared with 
the internal one. For these reasons, there is no direct clas- 
sical counterpart of an ultrarelativistic fluid and a relativistic 
description is needed even for small velocities. 

We solve the equations of relativistic hydrodynamics in a 
3D periodic domain usin g the high-resolution shoc k captur- 
ing scheme described in jRadice & Rezzollall2012l). In par - 
ticular, ours is a flux- vector-splitting sche me (Torolll999t). 
using the fifth-order MP5 reconstructi on dSuresh & Huvnhl 
1997), in local characteristic variables dHawkd 120011) . with 
a li nearized flux-split algorit hm with entropy and carbuncle 
fix jRadice & RezzollaKOll) . 

3. BASIC FLOW PROPERTIES. 

Our analysis is based on the study of four different mod- 
els, which we label as A, B, C and D, and which differ for the 
initial amplitude of the driving factor F = 1,2,5 for mod- 
els A-C, and F(t) = 10 + \t for the extreme model D. Each 
model was evolved using three different uniform resolutions 
of 128 3 , 256 3 and 512 3 grid zones over the same unit length- 
scale. As a result, model A is subsonic, model B is transonic 
and models C and D are instead supersonic. The spatial and 
time-averaged relativistic Mach numbers (vW) / (c s W s ) are 
0.362, 0.543, 1.003 and 1.759 for our models A, B, C and D, 
while the average Lorentz factors are 1.038, 1.085, 1.278 and 
1.732 respectively 

The initial conditions are simple: a constant energy density 
and a zero-velocity field. The forcing term, which is enabled 
at time t = 0, quickly accelerates the fluid, which becomes 
turbulent. By the time when we start to sample the data, i.e. at 
t = 10 (light-)crossing times, turbulence is fully developed 
and the flow has reached a stationary state. The evolution 
is then carried out up to time t = 40, thus providing data 
for 15, equally-spaced timeslices over 30 crossing times. As 
a representative indicator of the dynamics of the system, we 
show in the left panel of Fig.Q~]the time evolution of the aver- 
age Lorentz factor for the different models considered. Note 
that the Lorentz factor grows very rapidly during the first few 
crossing times and then settles to a quasi-stationary evolution. 



Furthermore, the average grows nonlinearly with the increase 
of the driving term, going from (W) ~ 1.04 for the subsonic 
model A, up to (W) ~ 1.73 for the most supersonic model D. 

Flow quantities such as the energy density, the Mach num- 
ber or the Lorentz factor show large spatial variance, even 
in our subsonic model. Similar deviations from the aver- 
age mass density, have been reported als o in classical turbu - 
lent flows of weakly compressible fluids dBenzi et al.ll2008l) . 
where it was noticed that compressible effects, leading to the 
formation of front-like structures in the density and entropy 
fields, cannot be neglected even at low Mach numbers. In the 
same way, relativistic effects in the kinematics of the fluid, 
suc h those due to nonlinear couplings via the Lorentz fac- 
tor dRezzolla & Zano tti 2002), have to be taken into account 
even when the average Lorentz factor is small. The proba- 
bility distribution functions (PDFs) of the Lorentz factor are 
shown in the right panel of Fig. Q]for the different models. 
Clearly, as the forcing is increased, the distribution widens, 
reaching Lorentz factors as large as W ~ 40 (i.e. to speeds 
v ~ 0.9997). Even in the most "classical" case A, the flow 
shows patches of fluid moving at ultrarelativistic speeds. Also 
shown in Fig. Q]is the logarithm of the Lorentz factor on the 
(y, z) plane and at t — 40 for model D, highlighting the large 
spatial variations of W and the formation of front-like struc- 
tures. 

4. UNIVERSALITY. 

As customary in studies of turbulence, we have analyzed 
the power spectrum of the velocity field 



E v {k) = \( \v{k)\ 
z J\k\=k 



dk . 



'|fe|=fc 

where fc is a wavenumber three-vector and 



v(k) 



v(x)e 



-2-irikx 



dx . 



(2) 



(3) 



with V being the three-volume of our computational domain. 
A number of recent studies have analyzed the scaling of the 
velocity power spectrum in the inertial range, that is, in the 
range in wavenumbers between the lengthscale of the prob- 
lem and th e scale at which d issipation dominates. More 
specifically, llnoue et all d201 lb has reported evidences of a 
Kolmogorov fc~ 5 / 3 scaling in a freely-decaying MHD turbu- 
lence, but has not provided a systematic convergence study 
of the spectrum. Evidences for a k~ 5 ^ 3 scaling were also 
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Figure 2. Power spectra of the velocity field. Different lines refer to the three 
resolutions used and to the different values of the driving force. The spectra 
are scaled assuming a fc -5 / 3 law. 

found bv lZhang et al. ( 2009), in the case of the kinetic-energy 
spectrum, which coincides with the vel ocity power-spectrum 
in the incompressible case. Finally, IZrake & M acFadyen 
(120121) has performed a significantly more systematic study 
for driven, transonic, MHD turbulence, but obtained only a 
very small (if any) coverage of the inertial range. 

The time-averaged velocity power spectra computed from 
our simulations are shown in Fig. [2] Different lines refer to 
the three different resolutions used, 128 3 (dash-dotted), 256 3 
(dashed) and 512 3 (solid lines), and to the different values of 
the driving force. To highlight the presence and extension of 
the inertial range, the spectra are scaled assuming a fc~ 5 ' 3 law, 
with curves at different resolutions shifted of a factor two or 
four, and nicely overlapping with the high-resolution one in 
the dissipation region. Overall, Fig. [2] convincingly demon- 
strates the good statistical convergence of our code and gives 
a strong support to the idea that the key pre diction of the Kol- 
mogorov model (K41) dKolmogorovl 199 lb carries over to the 
relativistic case. Indeed, not only does the velocity spectrum 
for our subsonic model A shows a region, of about a decade 
in length, where the fc~ 5 / 3 scaling holds, but this continues 
to be the case even as we increase the forcing and enter the 
regime of relativistic supersonic turbulence with model D. In 
this transition, the velocity spectrum in the inertial range, the 
range of lengthscales where the flow is scale-invariant, is sim- 
ply "shifted upwards" in a self-similar way, with a progressive 
flattening of the bottleneck region, the bump in the spectrum 
due to the non-linear dissipation introduced by our numerical 
scheme. Steeper scalings, such as the Burger one, k~ 2 , are 
also clearly incompatible with our data. 

All in all, this is one of our main results: the velocity power 
spectrum in the inertial range is universal, that is, insensitive 
to relativistic effects, at least in the subsonic and mildly su- 
personic cases. Note that this does not mean that relativistic 
effects are absent or can be neglected when modelling rela- 
tivistic turbulent flows. 

5. INTERMITTENCY. 

Not all of the information about relativistic turbulent flows 
is contained in the velocity power spectrum. Particularly im- 
portant in a relativistic context is the intermittency of the ve- 
locity field, that is, the local appearance of anomalous, short- 
lived flow features, which we have studied by looking at the 
parallel-structure functions of order p 

Sjjfr) = (\S r v\ p ), 5 r v=[v(x + r)-v(x)]-~ (4) 
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Figure 3. Compensated, third-order, parallel structure function computed for 
the different models as functions of r/A. Note the very good match with the 

classical s| ~ r behaviour. 

where r is a vector of length r and the average is over space 
and time. 

Figure [5] reports the compensated, third-order, parallel 

structure function, s\ , as functions of rj A, where A is the 
grid spacing. Within the inertial range, classical incompress- 
ible turbulence has a precise prediction: the Kolmogorov 4/5- 
law, for which ((S r v) 3 ) = |er, where e is the kinetic-energy 

dissipation rate. This translates into S3 ~ er. As shown in 
the figure, the structure functions are somewhat noisy at small 
scales, but are consistent with the classical prediction over a 
wide range of lengthscales, with linear fits showing deviations 
of ~ 5%, and an increase of e with the driving force. 

Although even in the classical compressible case, the 4/5- 
law is not strictly valid, we can use it to obtain a rough es- 
timate of the turbulent velocity dissipation rate dPorter et alJ 

120021) . We find that e, as measured from S 3 ' or directly from 
((8 r v) 3 }, grows linearly with the Lorentz factor, in contrast 
with the classical theory, where it is known to be independent 
of the Reynolds number. This is consistent with the observa- 
tions that in a relativistic regi me the turbulent ve l ocity shows 
an ex ponential decay in time (Zrake et afll201 lb llnoue etafl 
2011), as opposed to the power-law decay seen in classical 
compressible and incompressible turbulence. 
The scaling exponents of the parallel structure functions, 

(J have been computed up t o p = 10 using the extended-self- 
similarity (ESS) technique dBenzi et al.lll993l) and are sum- 
marized in Table [T] The errors are estimated by comput- 
ing the exponents without the ESS or using only the data at 
the final time. We also show the values as computed us- 
ing the classical K41 theo ry, as well as using th e estimates 
by She and Leveque (SL) (IShe & Levequd[T994l) for incom- 



pressible tur bulence, i.e. Qo 



2 - 2(|)p/ 3 , and those 



by Boldyrev (lBoldvrevl2002f) for Kolmogorov-Burgers super- 



sonic turbulence, i.e. Q 



+ l-(i)P/ 3 . 



Not surprisingly, as the flow becomes supersonic, the high- 
order exponents tend to flatten out and be compatible with the 
Boldyrev scaling, as the most singular velocity structures be- 
come two-dimensional shock waves. (J, instead, is compati- 
ble with the She-Leveque model even in the supersonic case. 
This is consistent with the observed scaling of the velocity 
power spectrum, which presents only small intermittency cor- 
rections to the A: -5 / 3 scaling. Previous classical studies of 
weakly compressible dBenzi et alJ 120081) and supersonic tur- 
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Table 1 

Scaling exponents of the parallel structure functions computed using the ESS technique and analytical predictions from the KS41, SL and Boldyrev models. 



Model 


cf 


^2 


C 11 

^3 


ci 


cl 


Cl 


Ct 


cl 


Cl 




K41 


0.33 


0.67 


1 


1.33 


1.67 


2 


2.33 


2.67 


3 


3.33 


SL 


0.36 


0.70 


1 


1.28 


1.54 


1.78 


2.00 


2.21 


2.41 


2.59 


Boldyrev 


0.41 


0.74 


1 


1.21 


1.39 


1.56 


1.70 


1.84 


1.96 


2.08 


A512 
B512 
C512 
D512 


0.37 ± 0.01 
0.36 ± 0.01 
0.37 ± 0.01 
0.38 ± 0.005 


0.70 ± 0.02 
0.70 ± 0.03 
0.70 ± 0.02 
0.71 ± 0.01 


1 ± 0.02 
1 ± 0.04 
1 ± 0.03 
1 ± 0.03 


1.27 ± 0.03 
1.27 ± 0.05 
1.26 ± 0.04 
1.25 ± 0.03 


1.51 ± 0.02 
1.50 ± 0.07 
1.48 ± 0.05 
1.46 ± 0.05 


1.72 ± 0.03 
1.70 ± 0.08 
1.68 ± 0.07 
1.64 ± 0.07 


1.89 ± 0.04 
1.86 ± 0.12 
1.84 ± 0.09 
1.79 ± 0.09 


2.04 ± 0.04 
1.99 ± 0.16 
1.98 ± 0.11 
1.92 ± 0.11 


2.17 ± 0.03 
2.10 ± 0.21 
2.09 ± 0.13 
2.04 ± 0.14 


2.27 ± 0.02 

2.18 ± 0.26 

2.19 ± 0.16 
2.14 ± 0.16 
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Figure 4. PDFs of the velocity v z for the different models considered (solid 
lines). As the forcing is increased, the PDFs flatten, while constrained to 
be in (—1, 1) (shaded area). Increasingly large deviations from Gaussianity 
(dashed lines) appear in the relativistic regime. 



bulence dPorter et alj|2002l) found the scaling exponents to be 
in very good agreement with the ones of the incompressible 
case and to be well described by the SL model. This is very 
different from what we observe even in our subsonic model A, 
in which the exponents are significantly flatter than in the SL 
model, suggesting a stronger intermittency correction. This 
deviation is another important result of our simulations. 

One non-classical source of intermittency is the genuinely 
relativistic constraint that the velocity field cannot be Gaus- 
sian as the PDFs must have compact support in ( — 1,1). This 
is shown by the behaviour of the PDFs of v z and plotted as 
solid lines in the shaded area of Fig. [4] Clearly, as the Lorentz 
factor increases, the PDFs become flatter and, as a conse- 
quence, the velocity field shows larger deviations from Gaus- 
sianity (dashed lines). Stated differently, relativistic turbu- 
lence is significantly more intermittent than its classical coun- 
terpart. 

6. CONCLUSIONS. 

Using a series of high-order direct numerical simulations 
of driven relativistic turbulence in a hot plasma, we have ex- 
plored the statistical properties of relativistic turbulent flows 



with average Mach numbers ranging from 0.4 to 1.7 and aver- 
age Lorentz factors up to 1.7. We have found that relativistic 
effects enhance significantly the intermittency of the flow and 
affect the high-order statistics of the velocity field. Neverthe- 
less, the low-order statistics appear to be universal, i.e. inde- 
pendent from the Lorentz factor, and in good agreement with 
the classical Kolmogorov theory. 
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